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A STATE  VARIABLE  APPROACH  TO  STATIONARY  TIME  SERIES 
WITH  APPLICATION  FOR  PREDICTING  DISPLACEMENT 
IN  A LARGE  UNDERWATER  SUSPENDED  ARRAY 


INTRODUCTION 


The  original  extensive  treatment  of  stationary  time  series  analysis  was 
developed  by  Wiener  (1949)  in  his  classical  work  "Extrapolation,  Interpola- 
tion and  Smoothing  of  Stationary  Time  Series,"  reference  1.  Wiener's  results 
were  expressed  in  the  frequency  domain  and  could  not  readily  be  extended  to 
nonstatior.arv  problems.  In  1960  Kalman,  reference  2,  extended  Wiener's  work 
to  nonstationary  time  domain  problems.  Although  the  real  power  of  the  Kalman 
filter  theory  lies  in  its  ability  to  handle  the  nonstationary  case,  computa- 
tional advantages  of  the  Kalman  filter  make  it  worth  considering  for  many 
stationary  problems. 


The  background  information  necessary  to  understand  the  Kalman  filter 
computational  techniques  has  been  included  in  this  report.  In  most  cases, 
equations  are  stated  without  proof,  since  lengthy  derivations  are  usually 
required  for  their  development.  Several  texts  (references  3,  4,  and  5), 
which  treat  this  subject  with  much  more  rigor  and  detail,  are  available. 
However,  most  of  the  development  and  notation  used  here  is  consistent  with 
that  used  by  Sage  and  Melsa,  reference  4. 


As  an  example  of  the  application  of  Kalman  filter  theory  to  prediction, 
the  displacement  fluctuations  of  a large  underwater  suspended  array  are 
investigated.  Both  the  1-  and  2-pole  prediction  estimates  along  with  the 
associated  error  variance,  are  determined. 


Prediction  of  the  array  displacement  is  motivated  by  a desire  both  to 
reduce  the  corruptive  effects  of  sampling  the  array  position  and  to  decrease 
the  number  of  active  transmissions  in  the  vicinity  of  the  array. 


BACKGROUND 


Some  of  the  notations  used  in  this  report  may  not  be  familiar  to  all  the 
individuals  interested  in  this  subject.  Therefore,  this  section  has  been 
included  as  an  aid  to  understanding  the  material  that  is  to  follow.  The 
information  is  admittedly  brief,  and  those  desiring  more  detail  should  refer 


to  Sage  and  Melsa,  reference  4. 


STATE  VARIABLE  NOTATION 


Linear  systems  can  be  represented  in  a vector  form  known  as  the  state 
vector  differential  equation, 
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x(t)  = F(t)  x(t)  + G(t)  w(t), 

where  the  state  vector  x(t)  is  related  to  the  system  dependent  variable  and 
its  derivatives. 

The  solution  of  the  state  equation  is  given  by 

x(t)  = <t(t,  tQ)  x(tQ)  + $(t,  x)  G(t)  w(t)  dx, 

o 

where  G(t,t0)  (called  the  state  transition  matrix)  is  given  by  the  solution 
of 

d [<D(t,  tQ)]  = $(t,  tQ)  = F(t)  $(t,  tQ), 

dt 

with  the  initial  condition  $(t0,  tQ)  = I. 

If,  in  addition,  F(t)  is  time  invariant,  the  state  transition  matrix  is 
given  by 

•(',  ‘ “o1  * “1F  * - * Vi  f""1' 

when  the  eigenvalues  of  F are  distinct,  the  <*£  are  given  by  solving  the 
following  sets  of  equations: 


. X.  (t-t  ) 

, v n-1  1 o 

ao  + aixi  + + Vl  X1  = e 


“0  + alX2  + + Vl  X2  = 6 


n-1  X2(t-to) 


• • 

• • 


a0  + “l  Xn  + 


. X (t-t  ) 
.n-1  n o 

+ Vl  Xn  - 6 


KALMAN  FILTER  THEORY 


Suppose  we  have  a system  whose  behavior  can  be  described  by  the  state 
equation 


x(t)  = F(t)  x (t)  G(t)  w(t) 


(message  model), 
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where  w(t)  is  a zero-mean  white  noise  process  with  covariance 

cov  {w(t),  w(x)}  * Vw(t,  T)  = ¥w(t)  6(t-x) 
and  the  initio,  mean  and  variance  of  the  state  x(t)  is  known 

VJx(o)  = E {x  (o)}  Vx(o)  = var  {x(o)}. 

Further  suppose  that  some  measurements  are  made  which  are  related  to  the 
state  x (t)  as 

z(t)  “ H(t)  x(t)  + v(t)  (observation  model), 

where  v(t)  is  zero-mean  white  noise  with 

cov  {v(  t ) , v(T)}  = Vv(t,  T)  = tv(t)  (S(t-T) 
and  w and  v are  uncorrelated: 

cov  {w  (t),  v ( T) } = Vw  v(t,  T)  = 0. 

Based  on  the  foregoing  model,  the  Kalman  filter  desires  to  make  an 
estimate  of  the  state  x(t).  In  addition,  the  estimator  must  meet  the 
following  requirements: 

(1)  The  form  of  the  estimator  is  linear, 

£(t)  = a ft)  + f A(t,  T)  z(T)  dT 

— — o o — 

(2)  The  bias  error  of  the  estimator  is  zero, 

E {x( t ) - x( t ) } = E {xe( t ) } = 0. 

(3)  The  trace  of  the  variance  of  the  error  is  a minimum, 

tr  (var  {x(t)  - x(t)  }}  minimum. 

The  estimator  which  can  be  developed  to  meet  the  foregoing  requirements 
is  the  "optimum  linear  minimum-error-variance  sequential  state  estimator," 
or  more  commonly,  the  "Kalman  filter."  A complete  summary  of  the  Kalman 
filter  algorithm  is  shown  in  table  1 and  diagrammed  in  figure  1. 

STATIOHARY  PROCESSES 

Many  problems  can  be  considered  stationary  or  almost  stationary.  This 
fact,  together  with  the  knowledge  that  the  Kalman  filter  has  many 
computational  advantages  over  the  Wiener  filter,  makes  it  worth  considering 
the  Kalman  filter  for  the  solution  of  stationary  problems. 


3 


TR  5893 


Table  1.  Continuous  Kalman  Filter  Algorithm 


Model s 


Algorithms 


Mess  age 


F i 1 ter 


xftl  = F( t 1 xft)  + G( t 1 wft)  i(t)  = F( t 1 3ft)  + K(t)  [z  ft)  - H( t ) xft)] 


Observati  on 


:f  t ) = H(t)  xft)  + v(  t ) 


Kft)  = V ft)  HT(t)  * *(t) 

P V 

X 

Error  Variance 


V (t)  = Ff  t)  V (t)  + V (t)  F (t) 
e e e 

x x x 

-V  ft)  HT(t)  * _1 ft)  H(t)  V ftl 
e v e 

x - x 


+ Gft)  ^ ft)  G ft) 
w 


Prior  Statistics 

E iw't)i  = E { v(t)t  =0  E ixfO)}  = U x f0) 

CO  V {wfO,  w(t  )}  = Yw  (t)  6ft-T) 

cov  {w  ft),  v(T)}  = cov  {x(0),  wft)}  = cov  fx(0),  _v(t)} 

var  {x(0)t  = Vx  '01 

cov  { v f t ) , v(t  )}  = 'Vv  (t)  6 f t-T  ) 


Initial  Conditions 


<'0)  = E { x(0l}  = U fO) 


V (0)  = var  {x  f 0)}  = V (0) 
e — x 

x — 


(1)  The  message  and  observation  models  are  time  invariant, 


x(t)  = Fx(t)  + G w(t) 
z(t)  = Hx(t)  + v(t). 

(2)  The  input  and  measurement  noise  are  at  least  wide  sense  stationary, 
Vw( t , T)  = R^(T)  = \ 6(t) 

vv(t,  x)  = Rv(T)  = \ fi(T). 

(3)  The  observation  interval  began  at  t = -°°. 

Since  x(t)  is  stationary,  Ve  (t)  =*  Ve  (0)  = Vg^  = constant 
and  the  error  variance  equation  "becomes  — 

0 = FV  tV  FT-HT!  "*HV  + G 'i'  GT, 
e e v e w 

XX  - X — 

where  the  Kalman  gain  is  given  by 

K = Vex  HT  V 

and  the  Kalman  filter  becomes 

x(  t)  = Fx(  t)  + K £z(  t)  - H x(  t)]  x(  -“)  = 0. 

The  main  difficulty  with  this  stationary  algorithm  is  in  the  solution  of 
the  error  variance  equation  known  as  the  matrix  Riccati  equation.  This 
difficulty  arises  because  the  equation  is  quadratic  and  many  solutions  are 
possible.  The  usual  method  of  solution  involves  solving  the  nonstationary 
error  variance  equation  until  the  solution  reaches  a steady  state  value. 
Because  the  solution  represents  variance  quantities,  the  solution  must  also 
be  positive  definite. 

The  stationary  estimation  problem  was  first  solved  by  Wiener,  and  is 
represented  by  the  following  diagram. 


v(t)  Noise 


Wiener  showed  that  the  optimum  linear  minimum  error  variance  filter 
WQ(s),  when  ^(t)  and  v(t)  are  uncorrelated,  is  given  by 
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where 

PR  corresponds  to  physically  realizable, 

R (s)  = message  spectral  density  = 3 R^(l)  = 3 (t  t T,  t),  and  A(s) 
is  given  by 

Rz(s)  = R^Cs)  + RyCs)  = A(s)  AT(-s). 

The  matrix  A(s)  is  such  that  the  rational  function  det  [A(s)]  has  all  its 
poles  and  zeros  in  the  left  half  of  the  s-plane. 

The  equivalence  of  the  Kalman  and  Wiener  filters  can  be  shown  by  letting 

2(t)  = H x(t). 

In  this  case,  it  can  be  shown  that  the  message  spectral  density  becomes 
Ry  (s)  = H(SI-F)-1  G^w  G(-SI-Ft)_1Ht, 
and  the  optimum  Wiener  filter  is  given  by 
WQ(s)  = H( SI-F  + KH)-1K. 

For  the  particular  case  of  scalar  observations,  the  spectral  factoriza- 
tion becomes 


Ry(S)  = [Ry(s)]  + [Ry(s)] 


where  the  + and  - represent  left-half  and  right-half  plane  poles  and  zeros, 
respectively. 

In  particular,  if  we  write  the  spectrum  as  a function  of  rational  poly- 
nomials 


K(s)]+  = 


do  + a,  s + . . . + a S 

1 m 

80  + Bj  S + ...  + Sn 


m < n, 


then  the  coefficient  matrices  in  the  Kalman  model  become 


0 1 0 
0 0 1 


0 0 . 
0 0 . 
P0-B1 


Based  on  measurements  to  time  tj,  there  are  three  classes  of  state 
estimation.  Consider  the  following  diagram,  depicting  one  component  of  the 
state  vector  x(t). 


If  we  define  x(t  | t^)  = E {x(t)  | , 

the  classifications  are: 

\ * (1)  Filtering  (t  = tj), 


x(t  | t^ ) = x(t^  | t1 ) = x( t1 ) , 

(2)  Smoothing  (t  < t^), 

(3)  Predicting  (t  > t^). 

The  Kalman  filter,  which  estimates  the  current  state,  has  been  covered  in 
the  previous  section.  The  smoothing  problem  is  of  interest  when  improved 
estimates  are  required  for  previous  values  based  on  the  collection  of  more 
data.  The  prediction  problem  arises  when  some  future  value  must  be  estimated 
from  present  measurements. 

The  smoothing  problem  would  amount  to  a fairly  large  extension  of  this 
material  and  will  not  be  covered  here.  The  prediction  problem,  however, 
amounts  to  a simple  extension  of  the  Kalman  filter  and  a treatment  of  this 
problem  follows. 

The  state  at  time  t is  given  by 
x(t)  * 0 ( t , t.)  x(t.)  + <I>(t,  t)  G(t)  w(t)  dx. 

i 1 t: 

Taking  the  conditional  mean  of  both  sides,  we  obtain 
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x(t  | tj)  * E j x(t)  | z/t^)}  = $(t,  t^)  x (tj) 


+ /J  *(t,  T)  G(X)  E jw(T)  | z(t  )}  dt. 
1 


Since  w(T)  is  a white  noise  process, 
E (wCO  | z^C  1 1 ) } * E {w(  T)  } = 0 
the  predictor  equation  becomes 

x(t  | t^)  • *(t,  t^)  x (tj) 


for  T tj , 


for  t > tj . 


Other  forms  of  the  predictor  equation  can  be  obtained,  depending  on  how  t 
varies  with  respect  to  tj.  However,  all  forms  are  the  same  for  stationary 
processes . 

The  prediction  error  is  given  by 

x ft  I t,)  = x(t)  - x(t  I t.), 

— e 1 — — l 


and  the  variance  of  the  prediction  error  is 

Vx  ^ ^ c i ^ = var  | *e  (t  ^ 1 1 ^ } ’ 

V (t  I t,)  = *(t,  t.)  V (t.)  *T(t,  t.) 
x 1 ’lxl  1 


+ J*  *(t,  T)  G(T)  V (T)  GT(T)  <tT(t,  T)  dT. 

*1  - 

The  Kalman  filter  can  be  extended  to  handle  other  classes  of  problems. 

For  instance,  the  problem  of  nonwhite  noise  can  be  handled  by  augmenting  the 

state  vector  with  additional  states.  Also,  the  Kalman  filter  can  be  used  in 

the  solution  of  nonlinear  estimation  problems,  provided  that  certain  justifi- 
able assumptions  can  be  made.  T)iese  extensions  of  the  Kalman  filter  will  not 
be  covered  here,  since  they  are  not  applicable  to  the  problems  under  consid- 
eration. 


EXPERIMENT 

The  preceding  sections  have  still  not  made  it  clear  as  to  how  the  Kalman 
filter  techniques  can  be  applied  to  a practical  problem.  The  author  thought 
that  the  best  way  to  demonstrate  these  techniques  would  be  through  the  use  of 
an  example.  The  example  selected  arose  from  some  measurements  made  on  a 
large  underwater  suspended  array,  reference  6.  The  end  result  of  the  signal 
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processing  will  be  an  estimate  of  the  variance  of  the  error  associated  with 
predicting  future  array  displacements  based  on  current  measured  displace- 
ments . 


DATA  ACQUISITION 

Three  transponder  assemblies,  figure  2,  were  installed  by  the  ALCOA  R/V 
SEA  PROBE  on  the  horizontal  portion  of  the  underwater  suspended  array.  Other 
transponders  were  installed  on  the  ocean  bottom  both  inline  and  normal  to  the 
array.  A diagram  of  the  resulting  experimental  configuration  is  shown  in 
figure  3. 

The  procedure  for  data  collection  was  to  interrogate  an  array  transponder 
(PI)  with  a shipboard  transponder  (PO).  The  array  transponder  would,  in 
turn,  interrogate  the  bottom  transponders  (P7  and  P8),  which  would  reactivate 
the  array  transponder.  From  the  sequence  of  arrivals  received  onboard  the 
ship,  it  was  possible  to  determine  the  slant  range  between  array  and  bottom 
transponders.  The  slant  range  could,  in  turn,  be  converted  to  Cartesian 
coordinates . 

The  array,  having  a large  mass,  was  assumed  not  to  move  during  the 
sampling  sequence,  which  lasted  approximately  10  seconds.  Further,  the 
travel  time  variability  due  to  ship  drift  was  found  to  be  predictable,  and 
thus  could  be  removed  from  the  computations.  The  array  position  was  sampled 
once  each  minute,  which  was  more  than  sufficient  to  prevent  aliasing. 


DATA  PROCESSING 

A sample  of  the  raw  data  is  shown  in  figure  4.  The  different  symbols 
represent  data  collected  from  different  transponder  units.  In  order  that  a 
complete  data  set  could  be  obtained,  the  raw  data  were  smoothed  to  eliminate 
bad  data  points  and  to  fill  in  missing  data  points.  The  smoothed  data  were 
then  sampled  at  6 minute  intervals  to  form  the  time  series  from  which  further 
processing  would  be  accomplished. 

A time  series  plot  cf  the  resulting  smoothed  data  is  shown  in  figure  5. 
The  length  of  the  sample  is  262  hours,  and  both  the  relative  displacement 
fluctuations  inline  with  (Ax)  and  orthogonal  to  (Ay)  the  array  are 
shown.  The  significant  low  frequency  periods  appearing  in  the  time  series 
are  motions  due  to  inertial  effects  of  ocean  currents  in  the  vicinity  of  the 
array.  This  is  perhaps  better  illustrated  through  the  use  of  the  autocorre- 
lation function,  figure  6,  where  the  principal  period  of  the  displacement 
fluctuations  is  22.5  hours.  This  period  corresponds  exactly  to  the  inertial 
period  at  the  latitude  of  the  array. 

The  spatial  characteristics  of  the  array  can  be  observed  from  the  cross- 
correlation function,  figure  7,  between  transponder  units  1 and  2 at  the  ends 
of  the  array.  Here  the  uniformity  of  motion  across  the  array  is  demonstrated 
by  the  high  value  of  the  crosscorrelation  coefficient  at  zero  lag.  The 
crosscorrelation  plot  also  shows  clear  evidence  of  the  inertial  period  in  the 
displacement  fluctuations. 
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Figure  5.  Displacement  Time  Series 
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Histograms  of  the  displacements  orthogonal  to  and  inline  with  the  array 
are  shown  in  figures  8 and  9,  respectively.  The  displacements  inline  with 
the  array  are  relatively  small,  and  will  therefore  not  be  considered  in 
further  analysis.  The  standard  deviation  of  the  orthogonal  fluctuations  is 
34.6  ft,  and  will  play  an  important  role  in  the  Kalman  filter  analysis  of 
these  data. 

Additional  information,  which  is  necessary  before  analysis  of  these  data 
can  proceed,  is  contained  in  the  displacement  fluctuation  power  spectrum, 
figure  10.  This  power  spectrum  was  obtained  as  a Fourier  transformation  of 
the  displacement  autocorrelation  function.  A peak  exists  near  the  inertial 
frequency,  but  is  not  prominent  for  two  reasons.  First,  the  peak  is  sup- 
pressed by  the  logarithmic  scale,  which  was  used  to  enhance  any  energy  at  the 
higher  frequencies.  Second,  the  spectral  resolution  at  the  inertial  fre- 
quency is  poor;  a longer  sample  record  would  yield  higher  resolution. 
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KALMAN  FILTER  ANALYSIS 


Before  the  Kalman  filter  can  be  applied  to  the  problem  under  considera- 
tion, two  issues  have  to  be  considered.  First,  since  the  version  of  the 
Kalman  filter  developed  in  this  report  was  a continuous  model,  the  discrete 
power  spectrum  will  have  to  be  converted  to  an  equivalent  continuous  power 
spectrum.  Second,  the  continuous  power  spectrum  will  have  to  be  approximated 
by  some  form  of  rational  function. 

The  details  of  converting  the  spectrum  are  shown  in  appendix  A.  The 
initial  spectrum  is  the  discrete  one-sided  spectrum  of  figure  10.  The  result 
of  the  spectral  conversion  is  shown  in  figure  11  and  represents  one  half  of  a 
two-sided  continuous  spectrum. 

Fitting  the  power  spectrum  with  a rational  function  involves  fitting  the 
power  spectrum  with  a function  having  the  form 


R 

y 


(s ,m,n) 


2 

a + a„s  + 
o 2 

b + b.s2  + 
o 2 


2m 


m < n. 


It  is  not  an  easy  task  to  fit  the  spectrum  with  functions  having  high 
orders.  One  method  is  described  by  Sanathanan  and  Koerner,  reference  7.  It 
was  decided  to  limit  the  form  of  approximations  in  order  to  keep  the  results 
simple  enough  to  obtain  closed-form  expressions.  The  approximations  were 
therefore  restricted  to  a class  referred  to  as  lowpass  Butterworth  spectral 
approximations.  The  lowpass  Butterworth  family  is  described  by  functions  of 
the  following  form 


Ry(s,  n) 
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1 


16 


TR  5893 


where  n is  the  number  of  poles  in  either  half  of  the  s-plane.  Only  the  1- 
and  2-pole  approximations  are  considered  here,  since  obtaining  closed-form 
expressions  for  higher  orders  is  time  consuming. 

The  details  of  the  lowpass  Butterworth  spectral  approximation  are  covered 
in  appendix  B.  Briefly,  for  a given  n,  the  two  parameters  to  estimate  are 
U)c  and  P.  The  two  conditions  to  be  satisfied  were  selected  to  be: 


(1)  Total  spectral  energy  = measured  data  variance, 

(2)  The  frequency  .approximation  error  is  minimized  at  the  inertial 
f requency . 

Meeting  these  considerations  determines  the  parameters  u)c  and  P.  The 
resulting  1-  and  2-pole  spectral  approximations  are  shown  on.  figure  11,  where 
they  are  superimposed  on  the  continuous  measured  spectrum. 


All  information  is  now  available  for  the  application  of  the  Kalman  filter 
algorithm.  The.  details  of  these  computations  are  carried  out  in  appendix  C 
(1  pole)  and  appendix  D (2  pole).  In  both  cases  (1  and  2 pole),  the  esti- 
mator and  the  variance  of  the  error  of  the  estimator  are  determined  for  the 
Kalman  filter  and  the  Kalman  predictor. 

In  the  example  under  consideration,  the  noise  v(t)  is  very  small.  Thus, 
the  variance  of  the  error  associated  with  the  Kalman  filter  is  nearly  zero 
and  almost  perfect  current  estimates  of  array  position  can  be  made.  Predict- 
ing the  position  of  the  array  at  some  time  in  the  future,  however,  is  subject 
to  substantial  error.  This  can  be  observed  in  figure  12,  where  the  standard 
deviation  of  the  error  is  plotted  as  a function  of  time  for  both  predictors. 
Note  that  the  prediction  error  starts  off  near  zero  for  the  filtered  estimate 
and  approaches  the  standard  deviation  of  the  measured  data  (34.6  ft)  as  time 
increases . 

If  we  arbitrarily  define  prediction  time  as 


I 


Prediction  time 


t 

P 


=>  Tv  (t 

e 

y 


1 

e 


the  prediction  times  for  the  1-pole  and  2-pole  predictions  are  16  and  109 
minutes,  respectively.  Thus  it  can  be  seen  that  the  prediction  time  is  very 
sensitive  to  the  spectral  approximation.  The  results  should  not  be  surpris- 
ing, however,  since  the  2-pole  spectrum  has  less  high  frequency  fluctuation 
than  the  1-pole  spectrum.  Since  the  slope  of  the  true  spectrum  lies  between 
the  slope  of  the  1-  and  2-pole  approximations,  it  is  likely  that  the  true 
prediction  time  is  bounded  by  the  two  prediction  time  estimates.  The  true 
prediction  time,  however,  can  be  approached  only  by  making  higher  order 
approximations  to  the  power  spectrum. 
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SUMMARY  AND  CONCLUSIONS 


The  application  of  the  techniques  used  in  this  report  might  seem 
unnecessarily  cumbersome  for  the  treatment  of  a problem  of  this  nature. 
However,  this  example  was  used  <uly  to  demonstrate  the  use  of  the  Kalman 
method  in  filtering  and  predicting  stationary  time  series  data.  In  actual 
practice,  a discrete  version  of  the  Kalman  filter  would  probably  be  used,  and 
numerical  results  (rather  than  closed-form  expressions)  would  be  desired. 

Used  in  this  manner,  the  Kalman  filter  is  simple  to  implement  on  a digital 
computer,  and  is  computationally  efficient. 

The  chief  problems  in  using  these  methods  with  stationary  time  series 
data  are  (1)  solving  the  matrix  Riccati  equation,  and  (2)  fitting  the 
power  spectrum  with  a rational  function  of  high  enough  order  to  be  a true 
representation  of  the  real  spectrum. 

In  the  particular  example  used  here,  the  prediction  of  array  position  was 
shown  to  be  highly  sensitive  to  the  spectral  approximation.  Even  in  the  best 
case  considered  (2-pole  approximation),  the  prediction  time  (109  minutes)  was 
relatively  short.  Higher-order  spectral  approximations  would  improve  the 
accuracy  or  the  prediction  time  estimate  but  would  not  increase  the  predic- 
tion time.  Thus  the  possibility  of  predicting  array  positions  for  a long 
time  in  the  future  based  on  past  measurements  seems  doubtful,  unless  some 
deterministic  behavior  can  be  established. 
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Figure  8.  Displacement  Histogram  (Ay) 


Figure  9.  Displacement  Histogram  (Ax) 
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APPENDIX  A 

CONVERSION  OF  THE  DISCRETE  POWER  SPECTRUM 
TO  AN  EQUIVALENT  CONTINUOUS  POWER  SPECTRUM 


The  theory  developed  in  this  report  is  based  on  a continuous  filtering 
model.  Thus  it  was  necessary  to  convert  the  discrete  power  spectrum  to  a 
continuous  power  spectrum  in  order  to  apply  the  results. 

The  length  of  the  sample  was 

L » 262  hr. 

The  sampling  rate  was 

fs  = 10  samples/hr, 

which  corresponds  to  a sampling  interval  of 


= 0.1  hr. 


The  total  number  of  samples  are 
N = = 2620  samples. 

The  power  spectrum  shown  in  figure  10  was  obtained  by  taking  the  Fourier 
transform  of  the  autocorrelation  function  shown  in  figure  6.  The  number  of 
correlation  lag  values  was  selected  to  be 

m = 0.1  N, 

which  would  make  the  standard  error  of  the  spectral  estimates 


e -ypjP  = \Jo7l. 


The  length  of  the  two-sided  autocorrelation  function  was  then 


2(0.1  N) 


52.4  hr. 


A discrete  Fourier  transform  of  a record  of  length  T results  in  an  equivalent 
bandwidth  of 


Af  = £ 
T 


= 0.0191  cycle/hr. 


Each  value  plotted  in  figure  10  represents  10  times  the  log  of  the  one- 
sided spectral  estimate,  i.e., 
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Pk  = 10  log  X(k). 

The  sum  of  all  spectral  estimates  is  equal  to  the  variance  of  the  data. 

02  = T?62  = x(k>  = C 34 . 6 ) 2 
Z k=l 

with  fk  = k Af* 

To  convert  the  one-sided  discrete  spectrum  to  a two-sided  equivalent 
continuous  spectrum,  a bandwidth  correction  must  be  made: 


10  log  [p]  = 10  log 


[r( 


1 (Hz)  \ 

Af  ChFT/' 


The  continuous  spectrum  which  results  is  that  of  figure  11,  where  the 
spectral  units  are  dh/ f t2//Hz/2. 


A-2 
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APPENDIX  B 


LOWPASS  BUTTERWORTH  SPECTRAL  APPROXIMATION 


Assume  that  the  power  spectrum  of  the  message  process  can  be  approximated 
by  a rational  function  of  the  form 


P P 

R (s,n)  = — or  R (^,n)  = — 

v 2n  v , / 2n 

'-(H  ‘t-) 

v c ' c 


ST  - -p 


r 

o +Wc  “ 


Functions  of  this  type  are  known  as  the  Butterworth  family,  where  '^c  is 
the  3 dB  down  point  and  n is  the  number  of  poles  in  either  half  plane  of  the 
approximation. 

The  message  spectrum  and  the  measured  spectrum  are  related 


z(t)  and  z(t)  = y(t)  + v(t), 


where  y(t)  and  v(t)  are  uncorrelated,  and  v(t)  is  assumed  to  be  white  noise 
with 

cov  {v(t),  v(t) } = H'v  6(t-l). 

Then  for  all  s or  U), 

Rz(s)  = Ry( s ) + Vfv  or  Rz(o>)  = Ry(w)  + ^v. 

If  we  integrate  both  sides  over  the  frequency  range  for  which  the 
measured  power  spectrum  was  computed,  we  have 


R-l 
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2TTf 


CO  = 

m 


s (10/3600) 

= 2tt  — 


tt  rad 
360  s 


J+a>m  R (u>)  du  = o2  = -5-  / 

J ■»  z 2tt  j — co 


+0) 


+00 


-U)  z 
m 


1 + 


(?) 


2n 


1 

da)  + J ”*  Y dco 

2tt  j —co  v 

m 


Since  co,,,  » coc,  we  can  write 

Poo 


2 - 


2n  sin  (it/2n) 


+ T 


V IT 


The  above  expression  relates  P and  u>c,  since  all  other  parameters  are 
known. 

The  other  criterion  used  to  fit  the  Butterworth  spectrum  to  the  measured 
power  spectrum  was  to  minimize  the  error  of  the  spectra  at  the  peak  value  of 
the  measured  spectrum.  Since  we  are  using  only  one  value,  this  is  the  same 
as  minimizing  the  squared  error. 

Thus  with 
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To  minimize  the  error  at  U)  = u)p, 


Solving  for  ujc  in  terms  of  wp  results  in  the  following  value  for  uc: 
. . l/2n 

uj  = (2n-l)  a)  . 
c p 


Evaluation  of  the  spectral  approximation  for  the  1-  and  2-  pole  Butter- 
worth  spectra  results  in  the  following  parameter  values. 


I 


B-3/B-A 
Reverse  Blank 
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KALMAN  1-POLE  FILTER 


APPENDIX  C 


1-POLE  COMPUTATIONS 


xf  1 1 = Fxft)  + Gw(t) 
z(  t ) = Hx(t)  + v(t) 


Message  model 


O'lsprvati  op  model 


y(t)  = Hx(t) 


The  power  spectrum  model  is 


R (si  = 
v 


= [ys)]+ 


The  spectral  factorization  yields 


R (s) 
v 


[pl/?«  l * rp'/2»  i 

c cl 

<0  + S '<)  - s 

L c j Lc  J 


which  results  in  the  following  scalar  quantities  for  the  matrices: 


F = -u) 


T 1 /2 
H = ? ' w 


C = 1 


cov  { w(  t ) . wf  T ) } = ¥ 6(  t-T ) = 1 6 f t— T ) 

w 

cov  { v( t ) , v(t)}  = Yv  6 (t-T). 

The  filter  equation  is 


x(t)  = [F-KH]  x(t)  + Kz(t)  * Ax( t ) + Kz(t), 
which  has  a solution  given  by 


x(t)  ■ ♦ ft,—)  x (-<»)  + J1  <!>.  ft,  t)  K(t)  z(t)  dx, 


where 


x (-<*>)  = 0. 


C-l 
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L 


The  Kalman  gain  is  given  by 


K = V (o)  HT  f"1, 
e v 

x 


and  the  error  variance  is  given  by  the  solution  of  the  matrix  Riccati 
equati on 

,T 


0 = FV  (o)  + V (o)  FT  - V (o)  HT  ¥ 1 HV  (o)  + G V G , 

o m Pvt:  w 


where  V (o)  = V ft)  for  stationary  processes, 
e e 

x x 

Solving  the  equation,  we  find 


and 


K = 


.1/2 


ir  p ' 

1/2  ) 

IrM 

i 

Ir  p 

1/2 

-1 

The  solution  of  the  transition  matrix  is  given  by 


V'- T'  ■ *A<t'1)  ■ V ■ V 


where 


a = e 

o 


X(t-T) 


and  X is  the  eigenvalue  of  A given  by  solving 

| Xl  - A | =0, 
where 

1/2 


A = F-KH 


P 

- U)  1 + 

V 


Again  by  solving  the  equations,  we  find  that 


X = - u> 


-0) 


= e 


1/2 


(t-T) 


and  the  resulting  transition  matrix  is 
C-2 
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' 


y(t|tj)  = y(tx)  *p(t,  tp. 

The  solution  of  the  transition  matrix  is  given  by 
F(t-t,) 


*F  (t,  tj)  = e 1 - V-  V 


where 


_ X(t-t  ) 
a = e 1 
o 


and  X is  the  eigenvalue  of  F given  by  solving 

| Xl  - F | = 0 


Solving  the  equations,  we  find 


F = U)  , O = e 
c o 


-"c  “'h’ 


and  the  resulting  transition  matrix  is 


$F(t’  tl)  = 6 


-u  (t-t.) 
c 1 


Substituting  into  the  predictor  equation,  we  find  the  predicted  displacement 
is  given  by 


-0)  (t-t.) 

c 1 


for  t >_  t^  • 


y(t|tl)  = y(tj) 

The  variance  of  the  error  of  the  predicted  estimate  is  given  by  solving 


Vg  (t|tj)  = HVe  (t|tj)  HT, 


where 


Ve  (t,tl)  = *F(t»  tl)  Ve  (tl)  *F  (t’  tl) 


♦ *_  (t,  T)  G(T)  y (T)  GT(T)  $T(t,  T)  dT. 

tj  F w 


In  the  problem  under  consideration,  this  reduces  to 


-2uj  (t-t  ) . 

Ve  <t1*l)"  Ve  (tl}  e * “53 


1 - e 


-2u)  (t-t.) 
c 1 


C-4 


•1 


! 
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where 

x(-°°)  = 0. 

The  Kalman  gain  is  given  by 

K = Ve  (o)  HT  ^HV  (o)  , 
x — 

and  the  error  variance  is  given  by  the  solution  of  the  matrix  Riccati 
equation 

0 = FV  (o)  + V (o)  FT  - V (o)  HT  'T1  HV  (o)  + G'V  GT, 
e e e ve  w 

XX  X X — 


where  V (t)  = V (o)  for  stationary  processes, 
e e 

x x 


De  fine 


A 
S = 


' 

r . p 

1/4-l 

i + — 

V 

Solving  the  equations,  we  find 
S 


V (o)  = 

e 


p“: 


T 


VlV 


pW2 

L c 


s (1  + S + S2) 


and 


K = 


V2 


1 / 2, . 
p w 


1/2 


The  solution  of  the  transition  matrix  is  given  by 
4>A(t,  T)  = e A(t-T)  = aQi  + ctjA 
when  the  0l£  are  the  solution  of  the  set  of  equations 


“o  * °1  xi 


XL(t-T) 


D-2 
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Mt-T) 

a0  + “l  X2  = 6 ’ 

and  the  Xj^  are  the  eigenvalues  of  A which  we  obtain  by  solving 

I XI  - a|  =0, 
where 


A = F - KH 


72oj  S . 

v c 1 

2 2 „2  -Jh* 

-0)  -a)  S v c 


c c 


Solving  the  equations,  we  find 
-/2u)c  (1  + S) 


X1  " 


Cl  + j)  and  X^  = 


-fa  (1  + S) 


(1  - j), 


also 


-w  (1  + S)  (t  - t) 
c 


a = e 
0 


a 


+ sin 


J2 

1 ui  (I  + S) 
c 


Eu>  (1  + S 


a)  f 1 + S)  (t  - T) 
c 


] 


S)  (t  - T) 


—to  (1  + S)  < t - T) 
c 


[to  Cl  + S)  ft  - Til 

1* j 


and  the  resulting  transition  matrix  is 

V 


$A(t,  t)  = 


L A 


where 


(t,  T) 

*A 

‘11 

A12 

(t,  x) 

$A 

2 1 

A22 

~u)  (1  + S) 

Ct  - t) 

C 

<J>  C t,  T>  = e 
A11 


CO  ( 1 ♦ S)  (t  - t) 

COS 

c 

L / J 

'M)  [~ 


S)  (t  - T) 

T5 


D-3 
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V!1’ T)  ~ u>  ,<f*  s) 

12  c 


-a)  (1  + S)  (t  - T) 

c 


CO  ( 1 + 
c 

1 

V 


S)  (t  - T) 


“CO  (1  + S)  (t  - T) 

-72(0  (1  + sz)  [u,  (1  + S) 

♦.  ft,  T)  - ‘ e P ,i„ U 

*21  1 * s L fi 


S)  (t  - 1 


-to  (1  + S)  (t  - t) 
c 


f.  (t,  T)  = e 
A22 


CO  (1  + S)  (t  - t) 
c 


Substituting  the  values  determined  into  the  solution  for  the  state  equation, 
we  find  that  the  state  estimate  is 


A ~ 1 /o  o ~ 

Since  y(t)  = Hx(t)  * p ' co  x,(t), 
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the  estimated  displacement  is 


—to  ( 1 + S ) ( t - T ) 
c 


y(t ) =/2  u>  S e Jl 


cos 


0)  (1  + S)  (t 

c 


- T) 


+ sin 


0)  (1  + S)  (t  - t) 

c 


z(  T ) dx 


„ A I , P 

S = J 1 + 


tl/4 


-1 


Also,  since 

V (t)  = V (o)  = HV  (o)  HT  = P(J4  V (o)  , 

e e e c e 

y y x xu 

the  variance  of  the  error  in  the  estimated  position  is  given  by 


V ( t ) = V (o)  =/2  U)  Y 


C V 


p 

1/4 

. r 

1 ♦ y 

- 1 

V 

KALMAN  2-POLE  PREDICTION 

The  predicted  displacement  is  given  by 

y ( 1 1 tj)  = H x (tj  tj), 

where 

S (t|t  ) = $F  (t,  tj)  $ (tj). 

Since 


x (t)  = 


x j ( t )' 


>-*2 


(t) 


'9(t>/  pl/V' 

lift)/  P1/2a£ 


we  can  write  the  predictor  equation  as 

y ( 1 1 1 , ) = y ( 1 1 ) *n(t,  t,)  + y(  t , ) 4>,,(t,  t,). 
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The  solution  of  the  transition  matrix  is  given  by 

FCt-tj) 

$F(t>  ti^  = e = Oq1  + alF» 

where  the  ai  are  the  solution  of  the  set  of  equations 


a0  + alXl 


XjCt-tj) 


a0  + alX2  = 6 

and  the  Xj  are  the  eigenvalues  of  F which  are  obtained  by  solving 
| XI  ~ F | = 0. 

Solving  the  equations,  we  find  that 


~U)  -U) 

XL  = -jf-  [l  + j]  and  X2  = tj-  [l  - j]  5 

al  so, 
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Substituting  into  the  predictor  equation,  we  find  the  predicted  value  is 
given  by 


y ( 1 1 1 x ) = e 


-U) 

~F 

v 2 \ 


——  (t-t  ) + sin  (t-t.)  y 

A A 


'/T  u I . 

+ - — sin  (t_t.  ) y for  t > t . 

L c A J I " 

The  variance  of  the  error  of  the  predicted  estimate  is  given  by  solving 


Ve  (t| t1)  = HVe  (tltj)  H , 


where 


Ve  (t,tl)  = *F(t»  tl)  Ve  (tl)  *F  (t’  V 


+ J"!  Mt,  t)  G(T)  V (T)  GT(T)  <I)T(t,  X)  dT. 
t^  F w 

In  the  problem  under  consideration,  these  equations  simplify  to 


Ve  (tltl}  = ^c  Ve  (t|tl) 

y *n 


Ve  = V [Ve  (tl}  V.  (t’  tl)  + Ve  (tl}  V<C’  tl) 

X11  11  L X11  11  X12  12  J 

+ *F  (t’  tl)  [Ve  (tl)  $F  (t>  tl)  + Ve  (tl}  $F  (t’  tl) 

12  ex21  'll  ex22  *12  J 


+ I1  4>2  (t,  x)  dx. 

Jt\  F12 


K 
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